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GABOR FRAMES OF GAUSSIAN BEAMS FOR THE 
SCHRODINGER EQUATION 

MICHELE BERRA, MARTINA BULAI, ELENA CORDERO, AND FABIO NICOLA 


Abstract. The present paper is devoted to the semiclassical analysis of lin¬ 
ear Schrodinger equations from a Gabor frame perspective. We consider (time- 
dependent) smooth Hamiltonians with at most quadratic growth. Then we con¬ 
struct higher order parametrices for the corresponding Schrodinger equations by 
means of fi,-Gabor frames, as recently defined by M. de Gosson, and we pro¬ 
vide precise L 2 -estimates of their accuracy, in terms of the Planck constant K. 
Nonlinear parametrices, in the spirit of the nonlinear approximation, are also 
presented. Numerical experiments are exhibited to compare our results with the 
early literature. 


1. Introduction 


The goal of this paper is to construct asymptotic solutions for Schrodinger equa¬ 
tions 


( 1 ) 


ihdfU = H(t)u 
u( 0) = n 0 , 


by means of Gabor frames in the semi-classical regime (J% —>■ 0 + ). Here t e [0, T], the 

initial condition n 0 £ L 2 (R d ) and the quantum Hamiltonian H(t) is supposed to be 
the ft-Weyl quantization of the classical observable H(t, X ), with X = (x, £) G M 2d . 


1.1. Literature overview. There are many results about asymptotic solutions 
for partial differential equations (PDE’s), especially when the initial value is a 
wave packet, i.e. it is well localized in the physical space and it oscillates with an 
approximately constant frequency. In particular, if the initial profile is Gaussian 
(a coherent state), the solution will be highly concentrated along the classical 
trajectory, according to the correspondence principle. Such a semi-classical analysis 
for Schrodinger-type equations were widely studied in several papers, see e.g. [6j 
□3 H9l im EU m, 139] and the textbooks 13 EHEEJ EH®]. 
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The natural idea of this work is to decompose the initial value in Q by means 
of a ft-Gabor frame whose atoms are Gaussian coherent states, construct asymptotic 
solutions for each of them, a so-called Gaussian beam , and finally by superposition 
obtain the asymptotic solution to ([Tj). The main issues are the following: 

• Construction of a parametrix via Gabor frames 

• Estimates in L 2 for the parametrix and the error term 

• Numerical results. 

Despite the simplicity of the idea, we do not know a fully rigorous treatment 
of this matter. There are various attempts (see e.g. (2J and references therein) 
where however several arguments are carried out only at a heuristic level and with 
numerical experiments. The present paper is devoted to a rigorous study of these 
issues for a class of smooth Hamiltonians with at most quadratic growth and, unlike 
the previous work, we address from the beginning a finer analysis, that is higher 
order approximations: the approximate solution is searched as a (finite) sum of 
powers of h , and the order of approximation can be arbitrarily large. 

1.2. Notation and (h -)Gabor frames. To be explicit, let us fix some notation. 

The h-Weyl quantization of a function H on the phase space M. 2d is formally 
defined by 

(2) Hu{x) = Op™[H]u(x) = (2t Th)~ d [ e ih ~ 1{x - y)p H(°^-,p)u(y) dydp 

for every u in the Schwartz space S(R d ). The function H is called the h-Weyl 
symbol of H. For zq = (xo,Po) £ M 2d , we define by T h (zo) the Weyl operator 

( 3 ) f h (z 0 )u(x ) = Op^[e ih ~ 1{poX - Xop) }u{x) = e arl(j **-*° p °M u {x - Xq ). 

Such operator meets the definition of the so-called ft-Gabor frames, introduced in 
[19] as generalizations of Gabor frames. For a given lattice A in R 2d and a non-zero 
square integrable function ip (called window) on M. d the system 

g\p,A) = {T h {z)p : z G A} 

is called a /i-Gabor frame if it is a frame for L 2 (M d ), that is there exist constants 
0 < a < b such that 

(4) «ll/lli<ElTT*Wv)l 2 <f>ll/lll. v/ei 2 (r4). 

zS A 

In particular, when h = (27r) _1 , the operator T(z 0 ) := T*- 2 ^ ^Zo) is the so-called 
time-frequency (or phase-space) shift 

T(z 0 )f(x) = e -™P°*°e 27rwa 7(z - ^o) = T X0 M P0 f(x), z 0 = (x 0 ,p 0 ), 
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where translation and modulation operators are defined by 

T Xo f(x) = f(x-x o) and M p J(x) = e 2 ^ x f(x). 

In this case we recapture the standard definition of a Gabor frame (see the next 
section for more details). Since they first appearance in the fundamental paper by 
Duffin and Schaffer |20] on non-uniform sampling of band-limited functions, frames 
have been applied in many fields of mathematics and physics. In particular, Gabor 
frames have been widely used in signal analysis, time-frequency analysis, quantum 
physics. Recently Gabor frames have been successfully applied in the study of 
PDE’s. In [SI [TT] they have shown to provide optimally sparse representations for 
Schrodinger type propagators and in [f3] reveal to be an equally efficient tool for 
representing solutions to hyperbolic and parabolic-type differential equations with 
constant coefficients. More generally, wave packet analysis and almost diagonal- 
ization of pseudodifferential and Fourier integral operators by Gabor frames have 
been performed in [101 [El [El HSJ H3 |26l [35] . 

Pursuing the work on deformation of Gabor frames, that have been investigated 
by many authors Hi m ee m, gh, de Gosson in pjij uses the Schrodinger evolution 
to deform Gabor frames. 

This paper is intended to be in some sense complementary to the work of de 
Gosson [E], because we use the frames to construct an approximate propagator 
for the Schrodinger equation. We also adopt almost the same notation as in [III]. 

1.3. Main Results. Here we consider h-Gabor frames where the window function 
is the standard Gaussian 

(5) Mx) = Tr-^e-l*! 2 / 2 

and its rescaled version 

( 6 ) <j>*(x) = ( t Th)- d/4 e~ lxl2/ ^ h) . 

We define the coherent state centered at z G M. 2d the function 

(7) #(*) = f r ^)0 $(*). 

Consider the solution z t = ( x tl pt. ) to the Hamiltonian system 

(8) x t = d p H(t,x t ,p t ), fit = ~d x H(t,x t ,p t ) 
with initial value Zo = (xo,po) and define 

(9) HgHt,X)= A'€R m 

M=2 7 ' 

It is well-known that the solution to the corresponding operator Schrodinger equa¬ 
tion for h = 1 

( 10 ) 


id,S t (z„) = OpflHg>(t)]St(zo) So(zo) = I, 
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is a metaplectic operator S t (zo) corresponding to the symplectic matrix S t (zo) via 
the metaplectic representation mm- Following the works PI EH |39], a natural 
ansatz for asymptotic solutions to 0. modulo 0(h {N+1 '> /2 ), iV G M, where the 
initial value is the coherent state cut, that is 

1 , jihd t u = H(t)u 

1«(0) = C 

is provided by the Gaussian beam 


( 12 ) 



N 

e* s(t,z ° ) f h ( Zt )M h -i/ 2 S t ( Z o) ^ h j/ 2 bj(t, x)</>o(x). 


3=0 


Here the symmetrized action S is defined by 


(13) 


S(t,z o) 



with a being the standard symplectic form; the metaplectic operator M h - 1/2 is 
given by 

M h -i/ 2 f{x) = h- d/i f(h~ 1/2 x) 

and the functions bj(t,x ) are suitable polynomials in x with coefficients depending 
on t, zo (bo = 1), as we shall see in the sequel. 

The construction of the parametrix via Gabor frames, having the previous Gauss¬ 
ian beams as building blocks, is performed as follows. Set 


(14) h = 2tt h, 

and consider a h-Gabor frame Q n (4>o, h 1//2 A) with A = aZ d x /3Z d , a, (3 > 0. Let 
be a dual window in S(R d ) (see the next section for details). For iV > 0, t G [0, T], 
the parametrix to ([!]) is defined by 

(15) [{/«">(()/]((,■)= V 

zGhi/ZA 

Observe that U^(0)f = f. 

The following assumptions will be imposed throughout the paper. 


Assumption (H). Suppose that the symbol H(t,X ) is continuous with respect 
to (t, X) G [0,T] x R 2d and smooth in X, satisfying 

(16) \d%H(t,X)\ < C a , V|a| >2, X G R 2d , t G [0,T], 


This is our main result. 


GABOR FRAMES OF GAUSSIAN BEAMS FOR THE SCHRODINGER EQUATION 5 

Theorem 1.1. Under the Assumption (H) and with the above notation, there 
exists a constant C = C(T) such that, for every f G L 2 (R d ), 

(17) < GI/IUw Vt e (o,r] 

and 

(is) uihdt-mw^n 

L 2 (R d ) < Ch iN+3>/2 \\f\\ LHm Vi e [0, T], 

IfU(t ) denotes the exact propagator, for every f G L 2 (R d ), 

(19) ||(f/W(i)-f/(i))/|| iS(IIti) <Cifi( ,v +W 2 ||/|| 1S(Ka) Vi 6 [0,T], 

The pioneering papers in this spirit, for Hamiltonians H{t,x,p) = V(x) + p 2 / 2, 
come back to [30, IT1]. More general Hamiltonians were considered in [?> [39], which 
inspired this work. 

In the framework of nonlinear approximation we can also consider nonlinear 
parametrices, constructed as follows. 

Let rj > 0 be a threshold, and for / G L 2 (R d ) consider the index set 

A = {z 6 ft 1/2 A : |(/,fW>l > u. 

and the nonlinear operator 

d N, (*)[/](*.-)= v (/.f 1 w)#*(«,•)• 

In particular, for ■// = 0 we recover U^ N \t). 

In this case we attain the following issue. 

Theorem 1.2. Under the Assumption (H) and with the above notation, there 
exists a constant C = C(T) > 0 such that, for every f G L 2 (R d ), p > 0, 

( 20 ) l|d< A ' , (i)[/]||L W <c( £ K/,f*(zb' , >r) 1/2 vt 6(0,71 

and 

(2!) 

IKifts, - Sxt))t/‘ N) [/]IU ( K<) < cfi (,v+3)/2 ( V K/,f"(z) T , ‘>| 2 ) 1/2 V(6[0,T]. 

IfU(t ) denotes the exact propagator, for every f G L 2 (R d ), 

(22) ||f/< JV| (l)[/]-C/(()/|| 11( «, ) <c( V |(/,f"(z)7'f) 1/2 

+ C7^ +1 )/ 2 ( ^ |(/,| 2 ) 1/2 Vt G [0,T], 

zG A v ,f 
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A similar nonlinear parametrix was constructed in [36] in the case N = 0. In 
particular, the estimate (21) already appeared there (Theorem 5.1), but with an 
additional factor ^y\A v j\ in the righ-hand side: that estimate blows up when 7] —> 0 
if / has an infinite number of non zero Gabor coefficients (i.e. |A 0 j| = +oo), 

whereas we see that this is not the case in (21). In addition, the parametrix in 


was constructed by means of a truncated Gaussian, which introduces a further 
error in the estimate. 

As a byproduct of these techniques, in Section 4 we will also extend the weak 
deformation of frames result in m Proposition 18] to the case of higher order 
deformations. However, since the approach is perturbative in nature, our result 
just holds for th 1 / 2 small enough, and no longer for every t as in [T9] . 

We end up by recalling that Gaussian beam methods have been employed to 
obtain asymptotic solutions to hyperbolic PDE’s in [33] and hyperbolic systems in 

[21I31I331EZ102103]. 


1.4. Numerical Results. Finally, in Section [5] we provide some numerical exper¬ 
iments. We study the Cauchy problem Q for an Hamiltonian function H(t,x,p) 
of the form 

2 

(23) H(t,x,p) = V( x) + P ~, 

with an oscillating potential. This is a standard setting for the so-called generalized 
harmonic oscillator and it is perfectly suited to discuss the behavior of the method 
in the presence of a potential hill and a potential well. 

We develop numerical algorithms using MATLAB and the powerful LTFATj^] 
package, see [PUTT]. We address the long time propagation of the beams using the 
reinitialization procedure described in [36]. 


2. Preliminaries and time-frequency analysis tools 


We refer to [23] for an introduction to time-frequency concepts and in particular 
to tm for applications to Mathematical Physics. For sake of brevity, sometimes 
we write xy = x ■ y, the scalar product on and |t| 2 — t ■ t, for t G M rf . The 
brackets (f,g) denote the inner product of L 2 (R d ), i.e. (f,g) = f f(t)g(t)dt and 
||/1|| = (/,/). The Schwartz class is denoted by S(R d ). For 1 < p < oo, the space 
£ P (A) is the space of sequences a = {oaIasa c> n a lattice A, such that 


a Up 



< oo. 


Targe Time Frequency Analysis Toolbox http://ltfat.sourceforge.net/ 
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We write Sp(d, M) for the group of symplectic matrices on M 2d , i.e., A G Sp(d, M) 
if A is a 2d x 2d invertible matrix such that A T JA = J , where 


(24) 



The standard symplectic form on R 2d is 

(25) a(z,z')— (z'fJz, z,z'eR 2d . 


The metaplectic group is denoted by Mp(d). Consider S G Mp(d ) with covering 
projection ir h : S H> S G Sp(d, M). The appearance of the subscript h is due to the 
fact that to the h-dependent operator Vpf(x) = e ~ lPx ' x /( 2h ) f( x ) (chirp) corresponds 


the projection 7 r s ( Vp ) = Vp, with Vp = 


, P = P T , and to the Fourier 


Id Orf 

—P 0 d / 

transform Jf{x) = (2mK)~ d P f Rd e~ lxx I h f(x') dx 1 corresponds n h ( J) = J, defined 
in (24). For details see [19, Appendix A] and the books [TZ, 18]. In particular, for 
A > 0 we shall use the metaplectic operator M\ G Mp(d ) defined (up to a sign) by 

(26) M x f(x) = X d/2 f(Xx), f G L 2 (R d ) 

and whose projection is n h (M\) = M\ 1 the symplectic matrix 


(27) 


M x = 


x-'h 


0 d 

A Id- 


In the sequel we shall often use the fundamental symplectic covariance formula 
(28) f\z)S = ST^S-'z) S G Sp(d, M). 

2.1. h-Gabor Frames. The definition of a h-Gabor frame is already contained in 
the introduction. Consider a lattice A in R 2d . The Gabor system 

£(<A a ) = {P(z)v,z G A}, 

(recall that T{z) = 7 - ( 2?r ) J (^)) is a Gabor frame for L 2 (R d ) if there exist constants 
a, b > 0 such that for every / G L 2 (R d ) 


(29) 


<•11/111 <£KI>fUM 2 <&ll/ll 


z£ A 


If (29) holds, then there exists a 7 
Gi'J, A) is a frame for L 2 (l 


G L 2 (R d ) (so-called dual window), such that 
and every / G L 2 (R d ) can be expanded as 


1 = X)(/. GO'rfAb = D</. T(z)~i)T(z)y, 


( 30 ) 
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with unconditional convergence in L 2 (R d ). In particular, if the window <p is a 
Gaussian function, then there exists a dual window 7 that is smooth and well- 
localized, in particular 7 G S(R d ) (see [131, l2Tj). 

In what follows we investigate some useful properties which let us switch from a 
h-Gabor frame to a standard Gabor frame and vice-versa. To reach this goal, we 
define the dilation matrix D h and its inverse (T > ?1 ) _1 


(31) 


D h 


f Id 0 d\ 

\0d hl d . J 


(D h )~ l 


f Id 0 d \ 
\0 d h~'I d ) 


(recall that h = 2nh). 


Proposition 2.1. Let D h be the matrix defined in (31). The system G n (p, A) is a 
h-Gabor frame if and only ifG{p, (D n )~ 1 A) is a Gabor frame and the frame bounds 
are the same. Moreover, every dual window 7 of the Gabor frame G(p, (D h )~ l h) 
originates a h-Gabor frame (^( 7 , A), dual frame of G h {p, A). 


Proof. The first part is straightforward and follows the pattern of [19] Prop. 7]. 
Precisely, the system G r '{p, A) is a h-Gabor frame if and only if there exist positive 
constants a, b > 0 such that Q holds. Setting z = (x, 2rchp) G A if and only 
if (x,p) G (D h )~ l A, and using T h (x,2nhp) = T[x,p ), the inequalities Q are 
equivalent to 

«II/Il2< E I(/>tV)| 2 < *11/111, v/6L 2 (r“), 

z£(D h )~ 1 A 

as desired. 

Now, consider a dual window 7 G L 2 (R d ) of the Gabor frame G(p, ( D h )~ 1 A ). 
Then every / G L 2 (R d ) can be expanded as 

/= (f>T'(z)<p)f(z)'y 

z£(D h )~ 1 A 

{f,f h (x,2TThp)v)f H (x,2Trhp)~/ 

(x,p)£(D h )~ 1 A 

{x,p’)e a 

that is the system {T n (z), z G A} is a dual frame of the ft-Gabor frame G h (p, A). 

□ 


Given the Gaussian window fio defined in (J5]) and the lattice A = a7L d x fiTfi, 
with a, f3 > 0, de Gosson in pH, Prop. 12] shows that the system G h {f>o,h 1 ^ 2 A) 
is a fi-Gabor frame if and only if a(3 < 1, this means by the previous proposition 



GABOR FRAMES OF GAUSSIAN BEAMS FOR THE SCHRODINGER EQUATION 


9 


that the system 

(32) g($, ah 1/2 Z d x (3h~ 1/2 Z d ) 


is a Gabor frame if and only if af3 < 1. This frame will be used for the numerical 
experiments in Section [5j 

Using the same arguments as in the proof of Proposition |2.1 
following characterization: 


we obtain the 


Proposition 2.2. Let D n be the matrix defined in (31) and consider a Gabor 
frame Q(ip,A). The system Q( 7 , A) is a dual Gabor frame of G (if, A) if and only if 
G h ( r y,D h A) is a dual h-Gabor frame of the h-Gabor frame Q h (ip,D h A). Moreover, 
the frame bounds of Q(<p, A) and Q h (ip,D h A) are the same. 


We shall work with h-Gabor frames where both windows and lattices are rescaled. 
Their dual frames behave as follows. We set 

(33) f h (x ) = h~ d/4 (f(h~ 1/2 x) = M h -i/ 2 <f(x) = M {2nh) -i/ 2 ip(x). 

Proposition 2.3. Consider a Gabor frame Q(ip,A). The system Q( 7, A) is a dual 
Gabor frame of Q(ip, A) if and only if G h ( , y h , h 1 ^ 2 A) is a dual h- Gabor frame of 
the h-Gabor frame Q h (ip h , h 1 ^ 2 A). Moreover, the frame bounds of Q(ip, A) and 
G h (ip h , h 1 ! 2 A) are the same. 


Proof. Consider the metaplectic operator M h - 1/2 G Mp(d ) and its inverse M h 1/2 G 
Mp(d). For every / G L 2 (R d ), we have M h -i/ 2 M h i/ 2 f = f. Given the Gabor 
frame G(<p,A) with dual Gabor frame G( 7 , A), by Proposition 2.2 and using the 
boundedness of the metaplectic operators on L 2 (R d ) we can write, for every / G 
L 2 (R d ), 


f = Y. ow, 7 

z£D h A 

= M h - 1/2 (M h ii2f 1 T h (M hl/ 2z)ip}T h (M hl/ 2z) 1 

zeh 1 / 2 A 

= {f,M h - 1/a ? h (M hl/ 2z)<p)M h - i/ 2 f fi (M fc i / 2 z) 7 , 

zeh 1 / 2 A 

= Y U,T\z)f)f\z)i h , 

zeh 1 / 2 A 


where we used T h (z)M h - 1/2 = M h -i/ 2 T h (M h i /2 z), by the covariance formula (28) 
and the definition (33). This ends the proof of the equivalence of the dual frames. 
The proof of the frame bounds follows the argument of the first part of the proof 
of Proposition 2.1 □ 
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3. Bounds for the parametrix in L 2 

3.1. Preliminary remarks. We assume for the Hamiltonian H(t, X) the validity 

(O') 

of Assumption (H). In particular, we consider its second order Taylor term H~ 0 
at zo = (xo,po), as in (J9l) and the corresponding propagator S t (zo) in (flol). We can 

p-h 

also consider the operator S t (zo) dehned by 

(34) mS t \z 0 ) = H%\t)S t h (z 0 ) S 0 \z 0 ) = I. 

This operator is related to S t via the formula (omitting the dependence on Zo for 
simplicity) 

St = M h -i/2 St.Mfri/2. 

Indeed we have 

ihS t = ihd t M n -i/2S t M h i/2 

= hM n -i/20pi [H^ (t, x, p)]S t M h i /2 
= HOpf[H^ 2 \t, h~ 1/2 x , h l/2 p)]M h - 1 /2 S t M h i /2 
= Op ™[H (2 \t,x,hp)}S t h 
= H&Tt)S t h . 

Remark 3.1. The action of S t (zo) or S t (zo) on a mudulated Gaussian function 
can be written down explicitly, for details we refer to Section [5] and the references 
quoted there. 

Remark 3.2. The projection S t (z 0 ) represents the flow of the linear system with 
Hamiltonian HzJ{t,X). As a matrix, S t (z 0 ) is in Sp(d,M), for every t G [0,T]. 
The entries of the matrix S t (z 0 ) depend on t e [0, T\ and z 0 G M 2d but they are 
bounded, because this is true for the coefficients of the polynomial H^(t, A") by the 
assumption (16). The same is true for the entries of the inverse matrix Sf l (zo). 

3.2. Evolution of a coherent state. Consider the Cauchy problem ( [IT] ). Let z t 

be the trajectory of the corresponding Hamiltonian system, with initial condition 
Zq- We will show that an approximate solution (Gaussian beam) is given by 


(35) 


^i 0 (^-) = e^ ) f ft (^)S%o)^ 


= eV( t *> )f R (2 t )M fi - 1/2 5 t (2 0 )M Rl/2 ^ 
= e i<5(Mo )f*( Zt )M R - 1/2 5 t (z o )0 o , 


where 5(t,Zo ) is the symmetrized action dehned in (13). More generally, we will 
consider higher order approximations 4^{t ) in the form (12). 
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To this end we will have to estimate the remainder term 

(36) ■) := ( ihd, - H(t■). 

The following computation of were carried out in [,3Ih 21] for Hamiltonians 
of the form H(t,x,p ) = V{x) + p 2 /2 and in |7, 39] for more general Hamiltonians 
with polynomial growth. Here we briefly sketch the main points for the benefit of 
the reader, because the formula given in [39l (70)] for R.i^ contains a number of 
misprints. An explicit computation show that 

(37) ihdte** {t ’ zo) = z t ) - H{t, z,)) 

= ~e^ t,Zo) (^(p t x t -p t x t ) - H(t,z t )). 

On the other hand, given a function /(x) we have 
(3S) 

ihd t T h (z t ) f (x) = 2 ) p tX + P tXt + PtXt )f( x - Xt ) - ihx t Vf(x - x t )^j 

= ' \ z t) y-pt'X- Pt.xt H--- ihxt V J / (x) 

= f n (z t )( PtXt ~ PtXt +xd x H(t,z t ) + hd p H(t,z t )(-iV))f(x), 

which implies 


(39) ihd t T\z t )M h -_ 1/a 


= T h (z t )M h -i / 2 


ihd t + 


p t x t - p t .x t 


+ h 1/2 xd x H(t , z t ) + h 1/2 d p H(t , z t )(-iV) 


We also have the covariance formula 


H(t)T\z t ) = T n (z t )0 P y[H(t, x + x t , h{p + h-Ppt)} 


which implies 

(40) mt)f\zt)M h -„2 = f h iz t )M H -i/20p^[Hit, h 1/2 x + x t , h x ' 2 p + p t )}. 
Finally a Taylor expansion yields 


(41) H(t, h 1/2 x + x tl h 1,2 p + p t ) = H{t, z t ) + h 1/2 xd x H(t, z t ) + H 1/2 pd p H{t , z t ) 

N +2 

+ hH®{t,X) + rt /2 H%{t,X) + ^ +3 )/ 2 r(f +3) (f,X), 

1 =3 


X) = Pffimt, z,)X\ X e R m 
hH 7 ‘ 


where 
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and 

r% +3> (t,X)= 1 Y. f d xH(t, Zt + 0ft 1/2 J>f)X>(l - ff) N+2 de. 

' >' | 7 |=Ar + 3^ 0 

By the formulas we obtain 

(m - m)<p% N (t , o = ( m - »o P n #<?(<. x)i 

1V+2 N 

- v fi'/ 2 o P ”(ffW(i,A')] - fi (»+3)/ 2 o P rHy 3 >((, x)]))s,( z „) v ^.(t, •)* 

1=3 j =0 

X+2 

= en s ^T\z t )M h - 1/2 S t (z 0 ) ( ihd t - 0 h l / 2 0^[H^{t, S t (X))\ 

1=3 

N 

- fi (N+3 >/ 2 0pr[^ +3 )((, S,(A))]) v 

j =0 


Now we choose b 0 (t,X) = 1 and bj(t,X), j > 1 solutions to 

fi<9t^-(f,X)0o(x) = ^i+fc=j+2 S'i(X))](6 fc (t,-)0 O )(^) 

\6 J (t,X) = 0. 


Remark 3.3. Since 0o zs a Gaussian function and Op f[H^(t,St(X))] are differ¬ 
ential operators with polynomial coefficients depending on t, zq, we see that bj(t, X) 
is a polynomial in X, having coefficients depending on t, zq which are bounded, for 
the entries of the matrix S t , as functions of t e [0, T], z 0 € R 2d , are bounded, as 
well as the coefficients of the polynomial H^(t,X), by the assumption (16). 


With this choice of bj(t,X ) we finally obtain the desired formula for R^: 


(42) R^\t,-)^(ihd t -Jm)cf> h z f(t,-) 

= V fi( I +‘)/ 2 o p ”[H<;>(i,5 t (A))]6 t (i,.)4 

Z+fc>iV+3 

3<l<N+2 

0<k<N 

N 

+ Y R< w+3+ ‘) /2 0 P r[r<y 3 >((, S,(X))]b t (t, ■)*o). 

/c=0 


3.3. Bounds for the parametrix: proof of Theorems |1.1| , |1.2( From now on 

we work with the lattice A = alL d x (3Z d , a, f3 > 0. We shall need the preliminary 
estimate below. 
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Theorem 3.4. Let Q h ((pQ, h 1 ^ 2 A) be a h-Gabor frame. There exists a constant 
C = C(T ) such that, for every sequence {c z : z G h 1 ^ 2 A} G l 2 (h 11,2 A) we have 


. 1/2 

2 ^ vt g [o,T] 


(43) || ^ V 

z£l i 1 / 2 A zeh 1 / 2 A 

and, 

(44) 

|| X -)IU 2 (Mrf) < C?i (iV+3)/2 ( X l c 

zGh 1 / 2 A z&h 1 / 2 A 


1/2 


Vt g [o,T], 


Proof. Let us prove (43). We have 

II X C ^’ Ar ( t ’-)lll2 (Kd) = X CzCP{(t) n z ' N {t,-),(t) h f (t,-)>- 

zGh 1 / 2 A ZjZ'eh 1 / 2 A 

We will prove that 

(45) |(^(t,-),^(t,-))l<^-V), 

where G G ^ 1 (/i 1 / 2 A), with ^ -norm independent of A By the Cauchy-Schwarz 
inequality in (, 2 {h l ^ 2 A) and Young inequality we have 

1/2 


XI C ^(V 

zGh 1 / 2 A 


< 


d ) - v X i c 

zGh 1 / 2 A 

< ll^ll^(/i 1 / 2 A) X l C; 

zGh 1 / 2 A 


X I X 

z£h 1 / 2 A z'^h 1 l 2 A 
2 




1/2 


It remains to prove (45). If we write down explicitly the expression of the functions 
cpf N (t, ■), f z ) v (t, •), we see that it is sufficient to prove that 

(46) i(r"( Jl )W./»s,( 2 )9i.r“(z;)M,- 1/a s,( z ')g 2 )| < g(z - /), 

with G G £ 1 (h 1 < ,2 A) as above, when gi,g 2 are Schwartz functions. 

If we denote by I z , z > the left-hand side of (46) we have, with z t = (x t ,p t ), z' t = 

(x'vP't), 

U,zf = I (f h (z t - z , t )M h -i / 2 S t (z)gi, M h -i / 2 S t (z')g 2 )\ 

= \(f h (h~ 1/2 (x t -x' t ),h 1/2 (p t - p' t ))S t (z)g 1 ,S t (z')g 2 )\ 

= \{T\h- l ' 2 {z t - z’ t ))S t (z)g i, S t (z’)g 2 )\. 

Here T 1 stands for T h with h = 1. 

Now, metaplectic operators are bounded d>(M d ) —> £(M d ) [T7] and, as already 
observed in Remark 3.2, the entries of the matrix S t (z) are bounded as functions 
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of t G [0, T\ and 0 G M 2d . Hence the functions S t (z)gi and S t {z')g 2 are Schwartz 
functions with Schwartz seminorms bounded with respect to t, z, z'. 

Since the map (essentially the Short-time Fourier transform in Time-frequency 
Analysis) 

is continuous S(R d ) x S(R d ) —* 5(M M ) |23], there exists a constant C = C(T) such 
that 


4, Z '<c , (i + r 1 / 2 |^-4|)-( 2d+1 ) 

< C( 1 + h~ 1,2 \z - z '|)- (M+1) 

where in the last inequality we used the fact that the maps z —>■ z t has an inverse 
which is globally Lipschitz continuous. This easily follows from the assumption 
(16) (see e.g. [9]). Hence (45) is verified with G(z) = (1 + h~ 1 ^ 2 \z\)^^ 2d+1 \ and 


IIG ? IU 1 (*. 1 / 2 a) — + 1^1) ^ d+1 ^ —C'< Too, 


5eA 


and C is of course independent of h. 


C z 


Let us now prove (]44|). The proof is similar to that of (43). Indeed, by arguing 

\t,x ) in (42) we see first of all that 


as above and using the expression for Rl 

we gain a further factor h ( iV + 3 V 2 i Moreover we are again reduced to estimate 


terms of the same form as in (46), where now the functions g\ and g 2 depend on 
t G [0, T], z,z' G M. 2d . The point is that their seminorms in the Schwartz space 


remain bounded. This follows from the fact that in (42) the functions bj(t,-)(j )o 
have seminorms in the Schwartz space bounded with respect to t G [0, T] and 
Zq G 


2d , and the pseudodifferential operators which act on them have symbols 

i.e. bounded together with all their derivatives) with 


in the Hormander class Sqq 


seminorms uniformly bounded with respect to t , Zq (which in turns is a consequence 


of Remark 3.2). This concludes the proof of (44). 


Proof of Theorem E3 The bounds ( [T7| ) and ( [l8| ) follow at once from ( |43| and ( |44| ) 
with 

c,= (/,fW>- 


Indeed, we have 


and 


/= E lt; 

z&h 1 / 2 A 

( E ai 2 ) 1 / 2 <cii/ia, 

2G/l 1 / 2 A 
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for a constant C > 0 independent of h, by Proposition 2.3 


Finally we prove (19). This follows from Duhamel formula: if U(t, s ) is the exact 
propagator, with U(s,s ) = /, and U(t) = U(t, 0), we have 

u w (t)f - U(t)f = -j x f U(t, s ) (ihd s - hTs j) U^ N \s)f ds, 


and (19) follows from (18) and the the fact that U(t,s ) is a unitary operator in 


Proof of Theorem L2. The estimates (20), (21) follow from Theorem 3.4 
Concerning (22) we see that now Duhamel’s formula reads 


u™= Uf- >(0) [/] - (7(0)/ - J‘ U(t, s)(itti, - H(s)) C/f >( S )[/] ds, 

which gives the desired result, using (|2T| and the representation 


yr ) (°)[/]- c, (°)/= e </t*(j) i h )€ 


This concludes the proof. 


□ 


4. Higher order deformation of frames 


With the above notation, consider the function 

N 

0 o ,N {t,x) = e^ s(t ' 0) f h (z t )S t h ( 0 ) ^2 h :/2 bj(t, h~ 1/2 x)(f) q(x). 

3=0 

Here z t is the integral curve of the Hamiltonian system ([8]) with initial condition 
z 0 = 0. Also, the bfs are constructed as in Section 3.2 by considering as initial 
value the Gaussian 0 q, which is centered at Zq = 0. Let f t be the Hamiltonian flow 
defined by H(t,X). In particular z t = f t { 0). 

The following result extends m Proposition 18]; it reduces to that result for 
N = 0, at least for th 1 ^ 2 small enough. 


Theorem 4.1. Assume the validity of Assumption (H) and consider a h-Gabor 
frame Q h {<ff )1 h 1 ^ 2 A). Then there exists a constant e > 0, depending only o n H , 
the frame bounds of G h ((f>Q : h l ^ 2 A) (which are independent of h by Proposition 2.3) 
and A, such that for t G [0,T], th 1 / 2 < e, the system G h ((/>o N (t, ■), / t (/i 1//2 A)) is a 
h-Gabor frame, with frame bounds independent of h. 
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Proof. Let 


N 


0o N (t,x) = 1/2 ^)0o(^), 

3=0 

so that 

By arguing exactly as in PI Proposition 18], we see that it suffices to prove 
that G h (4>Q ,N (t, •), h 1//2 A) is a Gabor frame, at least for th 1//2 small enough. By 
assumption we have 

a \\f\\L 2 (R d ) — £ K/.r*w^)i 2 <6||/n 12M 

zeh 1 / 2 A 

for some 0 < a < b independent of h. 

Set 

N 

fjQ N {t,x) := fo ,N (t,x) - <j>l(x) = h~ 1/2 x)(i) n 0 (x) 

3 =1 

(we used the fact that b 0 (t,x) = 1). By the triangle inequality it is sufficient to 
prove that 

£ \{f,f\z)^ N (t,-)?< a -\\f\\i HKdy 

zeh 1 / 2 A 

Using the representation 

/= £ (l.f\z)^)f\z)<j>l 

zeh 1 / 2 A 

with being a dual window in 5(M d ), and using Young inequality in £ 2 (/i 1 / 2 A) we 
are reduced to prove that 

L,z' ■= \(T%z )4 f H (z f )^ N (t, -))| < G(z - z') 

for some G e £ 1 (h 1/2 A) with ||G||^ < a/(26). 

Now, setting 

N 

• 00 (t,x) := '^^h^ 2 bj{t 1 x)(j) 0 (x) = 
j =i 


and z = (t 0 ,Po); z ' — (^o,Po), we have 

G, Z ' = |(T n (^ - ^)0o,0o ,JV )l 


= h~ d / 2 


e l ° n ° x (j)o(h l / 2 (x — (x 0 — x' 0 ))) , Go(h~ 1 / 2 x) dx 


7 ^Z 2 ! - (rf. - rf.’ 


= |( T l {h l/2 {z - ^0)00) 0 o))|• 
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Observe that ipo is a Schwartz function whose seminorms are dominated by th 1//2 , 
because 6^(0, x) = 0 for j > 1. Hence 

I z ,z' < Cth 1/2 ( 1 + h~ l ' 2 \z - z '|)- (2d+1 ), 

for some constant C = C(T ) > 0. If we define G(z) = Cth 1 ^ 2 ( 1 + /r _1 / 2 |^|) _ ( M+1 \ 
we have 

Y G(z) = Y Cth~^ 2 { 1 + |5|)-( 2d+1 ) < C'th 1/2 . 

z£h 1 / 2 A z&A 

The desired conclusion then follows if we choose ( C'th 1 Z 2 ) 2 < a/ (26). □ 


5 . Numerical Results 


The aim of this section is to construct a parametrix of order N = 0 for an 


Hamiltonian function H(t,x,p ) of the form (23), where 

V e C°°(R d ) with \d*V(x)\ < C a , V|a| > 2 , x G R d , 
so that Assumption (H) (cf. (Jl6|) is satisfied. 


5.1. Time Evolution of Gaussian Beams. In what follows we first recall the 
set of ordinary differential equations controlling the time evolution of the Gaussian 
beam for the quadratic time dependent Hamiltonian (t,X) defined in ([9]). The 
details and the proofs of what follows can be found in many works, we refer for 
instance to 0 Chapter 3]. Let us represent the quadratic form Hzl\t,X), X = 
{x,p), in 0 as 

(t, X) = ^ (K t x ■ x + 2 L t x -p + Gtp-p) 
where K t ,L t ,Gt are real d x d matrices, G f , K t being symmetric. Setting 


m 


K t LJ\ 
L t G t )’ 


where Lj is the transpose matrix of L t , the classical motion driven by the Hamil- 
(‘ 2 ') 

tonian H Zq is given by the Hamilton equation 


(47) 



The classical flow S t (zo) = G Sp(d, M) for the Hamiltonian H~o\t,X ) 

(see (10)) satisfies 

£^ 0 ) = JF(t)S t (z 0 ), S 0 = I d . 
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We now focus on the ^-dependent equation (34) and consider the Gaussian beam 
defined in (35). We start with the ansatz 

(48) S t h (z 0 )<j)o(x ) = a(t)e^ FtX ' x , 

where T t G E f /, the Siegel space of complex symmetric matrices V such that Of > 0 
(for details we refer, e.g., to [22, Chapter 5]) and a{t ) is a complex valued time 
dependent function. Observe that W must satisfy a Riccati equation and a(t ) a 
linear differential equation. Indeed, imposing that the right-hand side ipt( x ) — 
a(t)e^ TtXX of (48) satishes the equation 

(49) ihd t ^ t (x) = Hj$(t)i/> t (x), ^ 0 (x) = $}(a;), 

it follows that the matrices T t fulfill the following Riccati equation 

(50) r t = -K-r t L t -Ljr t -r t Gr t , 

with the initial conditions Tq = il d , whereas the function a(t) satishes 


(51) 


a(t) = -Rr (if + G,r,)a(t) 


with initial condition a(0) = (nh) d//4 . Let us introduce the matrices M t = A t +iB t , 
N t — Ct + iDt which are nonsingular (see [22, Chapter 5]) and fulfill by MR) 

(52) M t = L t M t + G t N t 

(53) N t = -K t M t -LjN t: 

M 0 = I d , N 0 = il d . Furthermore, it can be proved the equality 


r t = N t M t 


-1 


The solution of (51) can be computed easily. Indeed, using (52), we observe that 

Tr (Lf + G t T t ) = Tr(M t M t _1 ) 


and the Liouvillc formula 


^log(detM,) = Tr(M t M t - 1 ), 


yields at once the solution 

(54) a(t) = (7rh) _a!/,4 (det Mi) _1//2 . 

Finally, we rephrase (J35| as 

(55) <t>%%x) = ef 5 ^e>( x -^(S t \z 0 )^)(x - x t ), 
where 


S{t,zo) = -x 0 po + 


p s x s - H(s, z s ) ds 
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We now come back to the example (23) for which 


m = 


(d 2 x V(x t ) 0 d 

I o d i d 


The equations (18]), (52), (53), (56) become 


(57) 


x t = Pt 

Pt = —d x V (x t ) 

Mt = N t 

Nt = -d 2 x V(x t )M t 
S(t,z 0 ) = - V{x t ). 


with the initial conditions 


{x 0l p 0l M 0 = I d ,N 0 = iI d ,5{0,Zo) = 


These equations characterize the Gaussian beams (55). 

5.2. Construction of the parametrix. We consider a h-Gabor frame G h ((f>o, h 1 ^ 2 A), 
with A = aL d x /3Z d , a,(3 > 0. Let 7 h be a dual window in <S(M d ). Using Theorem 
0 an approximate solution with N = 0 to the Cauchy problem 

( 5g ) f it&tu = {V(x) - |A )u 

\tt(0 )=Uo, 


is provided by the expansion ( [15] ) : 
(59) 


[(7 (0) (() t<o](i,i)= (u 0 ,f h (z)"/ h )<l>1’°(t,x). 

z&h 1 / 2 A 


In the framework of nonlinear approximation, dealt in Theorem 1.2 , we can fix a 
tolerance rj > 0 , and consider the set: 

(60) A VjUo = {z G h 1/2 A : \(u 0 ,f h (z)^ h )\ > r]}. 

In this case our parametrix becomes 

U^\t)[uo\(t,x) = c s <$?{t,x). 


zeA- 
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5.3. Algorithms. We work in dimension d — 1 and assume that the initial datum 
Uq is a signal of length L e 2N, defined on the periodized unit interval / = [0,1]. 
The space and frequency grids are defined as 

X = {0,1/ L, ...(L — 1 )/L}, Q = {—L/2, —{L — l)/2 ... 0 ... L/2 - 1}, 

(with periodic boundary conditions). For details and a complete exposition on the 
parameters involved in the algorithms we address to the LTFAT documentation in 
http: //It fat. sourceforge.net/. 


Algorithm 5.1 (Coefficients, via Discrete Gabor Transfrom). Consider a signal / 
of length L. 

1. Define a > 0 to be the length of the time shift. The space translation is 
a/L. 

2. Define M > 0 to be the number of channels. The frequency translation will 
be of length L/M. In order to have a frame, the density needs to be less 
than 1, hence the parameters a and M must be chosen such that aM < 1. 

3. Compute the dilated gaussian window g = (71 h)~ 1 / i e~ x (recall d = 

1). There is a LTFAT routine available, which gives a periodic normalized 
Gaussian: g = pgauss(L, Lhn). 

4. Calculate the dual window. The LTFAT command is gd =gabdual(g, a, M, L). 

5. Calculate the coefficients c of the signal / using the Discrete Gabor Trans¬ 
form c =dgt(/, ga, a, Ad). By default it gives the following: 


L—l 

c(m, n) = '52f(l)gd(l-an)e- 2 * Um ' M , 
1=0 


for n = 0,..., L/a — 1, m 
expansion 


0,..., M — 1. Since, by (55), we need the 


L—l 

c(m, n) = - an)e- 2 ^ l - an)m/M , 

1=0 


we apply the routine phaselock, i.e. c =phaselock(c). 


We can use this algorithm as initial step for our procedure. 

Algorithm 5.2 (Gaussian Beam, standard setting). Consider the initial condition 
Uq of length L. 

1. Use Algorithm 5T to compute the coefficients c of m 0 . 

2. Set a threshold r/ > 0, using thresh(c, rj). 

3. Set the initial values for the ODEs (cf.(|57|)). The initial displacement is 
Xo = an/L where, for symmetry reasons, 

n = -L/2 + 1,...,-1,0,1,..., L/2. 


(61) 
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This is not restrictive, since the ODE is solved modulus one. The initial 
momentum is po = 2'KhmL/M , with 

(62) m = -M/2 + 1,..., —1,0,1, M/2. 


Set 


h(0,zo) = ( anirhm)/M , M 0 = 1, N 0 = i. 


4. 

5. 


To obtain the normalization of pgauss we set a(0) = g(L). 

Solve the set of ODEs for each (m,n), where m ranges over the integers in 
(62) and n over the integers in (61). This can be done in MATLAB using 
a standard solver. In our implementation we used ode45. 

Construct the solution (f) [u 0 ] (t, x ) using the exp function of MATLAB. 


5.3.1. Large-time behavior. The unique feature of Gaussian Beams, or nearly coher¬ 
ent states, is that they are effective even in the presence of caustics. Nevertheless, 
in certain cases, the large-time behavior can be troublesome. Precisely, the smallest 
eigenvalues of the matrix T t can drop quickly and the corresponding Gaussian in 
(48) starts to spread. This leads to a drop of quality in our solution. This phenom¬ 
ena is studied in [3Bj in dimension d — 1. The authors relate the sign of the Hessian 


of the potential V (x) to the spreading of the Gaussian in time. When the Hessian 
is positive, i.e. the so-called potentil well, the matrix T* is bounded. When it 
is constant, T t shows a linear decay in time and when the Hessian is negative, T t 
decays exponentially. The latter case is named potential hill and it needs to be 
treated very carefully. Although this treatment is far from being conclusive for the 
general case, it suits our 1—dimensional problem. Hence, we follow their idea of 
reinitialization. 

We monitor the decay of T f and as soon as it drops under a certain tolerance, 
we stop the propagation at time t, say. We then compute the solution uj and use 
it as initial value for the evolution in the time interval [t, T]. Let us describe the 
algorithm in detail. 


Algorithm 5.3 (Gaussian Beam, reinitialization). Consider the initial condition 
Uq of lenght L. 

1. Set U r = u 0 , where u 0 is the initial value as before. 

2. Compute the Gabor coefficients of U r , using Algorithm |5.1[ 

3. Solve the ODEs, setting the initial values as in Algorithm [572} Use ‘Events’ 
inside odeset to monitor the decay of the fourth component. With this 
command, we can set a tolerance and, if the matrix T t drops below that, 
the computation of ode45 stops and return the time t together with the 
solutions at that time. 

4. Construct the solution U^ (t) [w 0 ] (t, x) at time t — t. 
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5. Set U r = uj) 0 \t){uo}(t, x) and execute again steps 1-4 until t reaches the 
final time. 


5.4. Numerical Results. The problems we consider are the ones presented in 
and, for sake of consistency, the same comparison method is used, i.e. the Strang 
Splitting pseudo-spectral method iHEj. 

5.4.1. Potential well. We consider the potential 

(63) V(x) = cos(27tx), 

treated in the Example 1 in [36j Subsec. 6.1.1] and set the initial values 

u 0 (x)=e- 25 ^-°^e 1 


iT p( x ) 

h 


(64) 


t„(x) =-- log(e (5(l -° 5 » + e - 5 (*-° 5 >) ; 
5 


with signal length L = 1024 and Plank constant h = -W-. We also set a threshold 


for the coefficients (cf. (60)) taking the ones with absolute value greater than 
r] = 0.01. With this tolerance, our initial reconstruction is still accurate. Indeed, 
at time t = 0 the relative error is just of order 0.3%, similar to the one in the 
Example 1 in PSl Subsec. 6.1.1]. 

The potential (63) is non negative for x G [0.25, 0.75]. So we expect the solution 
to be accurate inside this interval, since the beam width is constant. This is con¬ 
sistent with the results of Figure [lj where no reinitialization is needed. 


We can push our analysis a little bit further. In Figure [2] our initial values are 
represented by the picture on the left-hand side, whereas on the right-hand side 
we show the same function multiplied by a narrower Gaussian. This provides an 
“almost” compactly supported function inside the interval [0.25,0.75] where we 
have the potential well. In this case we obtain a 35% drop on the relative error. 
The solution is shown in Figure [3j which contains a comparison of the exact and 
beam solution. 


5.4.2. A potential Hill. Consider the initial value to be (64), we take the signal 
length to be L = 1024 and the Plank constant h = as above. We also 

set the same threshold 77 = 0.01. We consider the potential V(x) = cos(27 t(x + 
0.5)), whose Hessian is non-positive for x G [0.25,0.75]. This yields unacceptable 
numerical results for the standard algorithm, as shown in Figure [4| If we use the 
reinitialization algorithm, then our approximation improves a lot, see Figure [5j I 11 
this case we subdivide the time interval [0, 2] in eight uniform subintervals. We 
can try to do as before and chose an initial value which is compactly supported 
in [0,0.25] U [0.75,1], i.e. where the Hessian is non negative. The easiest way to 
do it is to shift the initial value used above by half the length of the signal, see 
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(c) Comparison, center. (d) Comparison, side. 

Figure 1. Potential V(x ) = cos(27nr). (A) The exact solution. (B) 
The beam solution. (C) and (D) Comparison of the exact and beam 
solution, “o”: the real solution; the beam solution. 




FIGURE 2. Initial values. The right one is numerically supported in 
[0.25,0.75] 
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(a) Strang-Splitting solution 


(b) Beam solution 




(c) Comparison, center. 


(d) Comparison, side. 


FIGURE 3. Potential V(x) = cos(27nr), compactly supported initial 
values. (A) The exact solution. (B) The beam solution. (C) and (D) 
Comparison of the exact and beam solution, “o”: the real solution; 
the beam solution. 


Figure [6j Then, as shown in Figure [7| even without any reinitialization we get a 
perfect reconstruction, as expected. We notice that this is nothing but the first 
case shifted by half the wavelength. 


5.4.3. A potential hill and well. Consider the initial value to be (64), we take the 
signal length to be L = 1024 and the Plank constant h = as above. We also 
set the same threshold rj = 0.01. 

We take V(x) = 10 + sin(27r(x + 0.5)), whose Hessian is non-positive for x € 
[0,0.5]. Once again, the numerical results without reinitialization are far from 
being consistent, see Figure [8} If we use the reinitialization, our results improve 
greatly as shown in Figure [9j 

We can try a similar trick as before and pick the usual compactly supported 
window shifted in the interval [0.5,1], once again the results are very good, see 


Figure 10 

















































GABOR FRAMES OF GAUSSIAN BEAMS FOR THE SCHRODINGER EQUATION 


25 



(a) Strang-Splitting solution 



(b) Beam solution 



(d) Comparison, side. 


Figure 4. Potential V(x) = cos(27r(a; + 0.5)). (A) The exact solu¬ 
tion. (B) The beam solution. (C) and (D) Comparison of the exact 
and beam solution, “o”: the real solution; the beam solution. 


5.5. Future Works. The LTFAT package provides a 2-dimensional discrete Gabor 
transform (dgt2), although the phaselock command is not yet available. Nev¬ 
ertheless, it seems possible to implement a multidimensional algorithm using the 
same approach. We plan to develop such a command and then investigate the two 
and three-dimensional cases. 
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(a) Strang-Splitting solution 



(c) Comparison, center. 



(b) Beam solution 



(d) Comparison, side. 


Figure 5. Potential V(x) = cos(27r(a; + 0.5)) with reinitialization. 
(A) The exact solution. (B) The beam solution. (C) and (D) Com¬ 
parison of the exact and beam solution, “o”: the real solution; 
the beam solution. 


References 

[1] G. Ascensi, H. G. Feichtinger, and N. Kaiblinger. Dilation of the Weyl symbol and Balian- 
Low theorem. Trans. Amer. Math. Soc., 366(7):3865-3880, 2014. 

[2] G. Bao, J. Lai, and J. Qian. Fast multiscale Gaussian beam methods for wave equations in 
bounded convex domains. J. Comput. Phys ., 261:36-64, 2014. 

[3] G. Bao, J. Qian, L. Ying, and H. Zhang. A convergent multiscale Gaussian-beam parametrix 
for the wave equation. Comm. Partial Differential Equations, 38(1):92 134, 2013. 

[4] W. Bao, S. Jin, and P. A. Markowich. On time-splitting spectral approximations for the 
Schrodinger equation in the semiclassical regime. J. Comput. Phys., 175(2):487-524, 2002. 

[5] W. Bao, S. Jin, and P. A. Markowich. Numerical study of time-splitting spectral discretiza¬ 
tions of nonlinear Schrodinger equations in the semiclassical regimes. SIAM J. Sci. Comput., 
25(l):27-64, 2003. 

[6] M. Combescure and D. Robert. Coherent states and applications in mathematical physics. 
Theoretical and Mathematical Physics. Springer, Dordrecht, 2012. 












































































































GABOR FRAMES OF GAUSSIAN BEAMS FOR THE SCHRODINGER EQUATION 


27 



Figure 6. Shifted initial value 



(a) Strang-Splitting solution 

Figure 7. Potential V ( x ) = cos 
initial values. 
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(c) Comparison, center. 


(d) Comparison, side. 
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exact and beam solution, “o”: the real solution; the beam 
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(c) Comparison, center. 


(d) Comparison, side. 
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the beam solution. 


[19] M. A. de Gosson. Hamiltonian deformations of gabor frames: First steps. Appl. Comput. 
Harmon. Anal., in press, 2014. 

[20] R. J. Duffin and A. C. Schaeffer. A class of nonharmonic fourier series. Trans. Amer. Math. 
Soc., 72:341 366, 1952. 

[21] H. G. Feichtinger and N. Kaiblinger. Varying the time-frequency lattice of Gabor frames. 
Trans. Amer. Math. Soc., 356(5):2001-2023 (electronic), 2004. 

[22] G. B. Folland. Harmonic analysis in phase space, volume 122 of Annals of Mathematics 
Studies. Princeton University Press, Princeton, NJ, 1989. 

[23] K. Grochenig. Foundations of time-frequency analysis. Applied and Numerical Harmonic 
Analysis. Birkhauser Boston, Inc., Boston, MA, 2001. 

[24] K. Grochenig and Y. Lyubarskii. Gabor frames with Hermite functions. C. R. Math. Acad. 
Sci. Paris, 344(3):157-162, 2007. 

[25] K. Grochenig, J. Ortega-Cerda, and J. L. Romero. Deformation of Gabor systems. ArXiv 
e-prints, Nov. 2013. 

[26] K. Grochenig and Z. Rzeszotnik. Banach algebras of pseudodifferential operators and their 
almost diagonalization. Ann. Inst. Fourier (Grenoble), 58(7):2279-2314, 2008. 
















































































30 


MICHELE BERRA, MARTINA BULAI, ELENA CORDERO, AND FABIO NICOLA 




(a) Strang-Splitting solution (b) Beam solution 

Figure 10. Potential V(x) = 10 + sin(27r(a; + 0.5)), compactly sup¬ 
ported initial values. 
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